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Abstract. The motto of this paper is: Let's face Bose-Einstein condensation through nonhnear 
dynamics. We do this by choosing variational forms of the condensate wave functions (of given 
symmetry classes), which convert the Bose-Einstein condensates via the time-dependent Gross- 
Pitaevskii equation into Hamiltonian systems that can be studied using the methods of nonlinear 
dynamics. We consider in particular cold quantum gases where long-range interactions between 
the neutral atoms are present, in addition to the conventional short-range contact interaction, viz. 
gravity-like interactions, and dipole-dipole interactions. The results obtained serve as a useful guide 
in the search for nonlinear dynamics effects in numerically exact quantum calculations for Bose- 
Einstein condensates. A main result is the prediction of the existence of stable islands as well as 
chaotic regions for excited states of dipolar condensates, which could be checked experimentally. 
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1. INTRODUCTION 

At sufficiently low temperatures a condensate of weakly interacting bosons can be rep- 
resented by a single wave function whose dynamics obeys the Gross-Pitaevskii equation 
[[I1I3. The Gross-Pitaevskii equation can be thought of as the Hartree equation for the 
ground state of interacting identical bosons, all occupying the same single-particle 
orbital ^f. Because of its nonlinearity the equation exhibits features not familiar from 
ordinary Schrodinger equations of quantum mechanics. For example, Huepe et al. 
demonstrated that for Bose-Einstein condensates with attractive contact interaction, de- 
scribed by a negative 5-wave scattering length, bifurcations of the stationary solutions 
of the Gross-Pitaevskii equation appear, and determined both the stable (elliptic) and 
the unstable (hyperbolic) branches of the solutions. The bifurcation points correspond 
to critical particle numbers, above which, for given strength of the attractive interaction, 
collapse of the condensate sets in. In Bose-Einstein condensates of ^Li [[51 [6l and ^^Rb 
atoms [iTliSI these collapses were experimentally observed. 

In those condensates the short-range contact interaction is the only interaction to 
be considered. In Bose-Einstein condensates of dipolar gases [(H HOl [HI [121 IHl also 
a long-range dipole-dipole interaction is present. Alternatively, following a proposal 
by O'Dell et al. [|14l [151 . by using a combination of 6 triads of appropriately tuned 
laser light condensates can be produced in which an attractive long-range gravity- 
like l/r interaction is present. These types of condensates offer the opportunity to 
study degenerate quantum gases with adjustable long-range and short-range interactions. 
While the experimental realization of condensates with gravity-like interaction lies still 
in the future, the achievement of Bose-Einstein condensation in a gas of chromium atoms 



[fT6]| . with a large dipole moment, has opened the way to promising experiments on 
dipolar gases [IVJ, which could show a wealth of novel phenomena fTE[ [T9l l20l |2T| . 
In particular, the experimental observation of the collapse of dipolar quantum gases has 
been reported [|22l| which occurs when the contact interaction is reduced, for a given 
particle number, below some critical value using a Feshbach resonance. 

In this experimental situation it is most timely and appropriate to extend the investi- 
gations of the nonlinearity effects of the Gross-Pitaevskii equation to quantum gases in 
which both the contact interaction and a long-range interaction is active, and this is the 
topic of the present paper. 



2. SCALING PROPERTIES OF THE GROSS-PITAEVSKII 
EQUATIONS WITH LONG-RANGE INTERACTIONS 

2.1. Gravity-like interaction, isotropic trap 

For an additional gravity-like long-range interaction Vu{r^r') = —u/\r — r'\ the time- 
dependent Gross-Pitaevskii equation for the orbital i// reads 

- A + fr^ +mn- 1 V/(r, t)\^-2N[ ^"^[^ '!}f d^T'] t) = i^ , (D 
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where we have used f23\ the "Bohr radius" a„ = h^/ (mu), the "Rydberg energy" Eu = 
fp'/{lma^), and the Rydberg time h/Eu as natural units of length, energy, and time, 
respectively, to bring the equation in dimensionless form. Furthermore, in ([T]) a is the 
5- wave scattering length, which characterizes the strength of the contact interaction V^. = 
5{f —r')AKafP' /m, N is the particle number, and y= ticOo/{2Eu) is the dimensionless 
trap frequency. It can be shown [23] that the solutions of ([T]) do not depend on all these 
three physical quantities but only on the two relevant parameters y/N^ and N^a / Uu- Thus 
one has, e.g., for the mean-field energy E{N,N^a* /au,y/N^)/N'^ = E{N = l,a/a„,7), 
with N^a* /uu = ajuu- 



2.2. Dipolar interaction, axisymmetric trap 



In Bose condensates of ^^Cr atoms |fT6l, which possess a large magnetic moment of 
jU = 6jUb, the long-range dipole-dipole interaction 
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must also be considered. Defining the dipole length by = jlQjl^m/ {InfP'), and using as 
unit of energy E^ = / (2ma^), of frequency cOd = 2E^/h and of time ti/E^, one obtains 
the Gross-Pitaevskii equation for dipolar gases in axisymmetric traps in dimensionless 



form 

= i^^\j/{r,t). (3) 

The physical parameters characterizing the condensates are the particle number A'^, the 
scattering length a /ad and the trap frequencies jp and perpendicular to and along the 

2/3 1/3 

direction of alignment of the dipoles (alternatively, the geometric mean (7 = 7p Yz and 
the aspect ratio A = y^/ Jp is used). However, a closer inspection of the scaling properties 
of ([b]) reveals [|24ll that the solutions depend only on three parameters, viz. A^^y, A, aja^- 
For the mean-field energy, e.g., the scaling law reads E{N ,a/ a^^N^T ■:^) = = 
X.ajad.y.X) /N, with A^^y* = y 



3. QUANTUM RESULTS: SOLUTIONS OF THE STATIONARY 
GROSS-PITAEVSKII EQUATIONS 

For the \/r interaction (monopolar quantum gases) we have solved [|23l the stationary 
Gross-Pitaevskii equation both variationally, using an isotropic Gaussian-type orbital 
1// = Aexp(— ^^r^/2), and numerically accurate, by outward integration of the nonlinear 
Schrodinger equation. For the dipole-dipole interaction (dipolar quantum gases) we have 
performed a variational calculation [|24| using an axisymmetric Gaussian-type orbital 
V/ = Aexp(-A:2p2/2 -fc,V/2). 

Fig.[T]shows the results for the chemical potential (eigenvalue of the stationary Gross- 
Pitaevskii equation) for the two interactions, plotted as a function of the scattering 
length. As can be seen, below a critical scattering length no stationary solutions exist, 
while two stationary solutions are born at the critical scattering length in a tangent 
bifurcation. At the bifurcation point the chemical potential, the mean-field energy, and 
the wave functions of the two branches of solution are identical. Such behavior is 
obviously a consequence of the nonlinearity of the underlying Schrodinger equation, 
and is reminiscent of exceptional points [|25ll26l discussed so far only in the context of 
open quantum systems with non-Hermitian Hamiltonians (see Ref. [27 J for references). 
In fact, a closer inspection shows [|27l that the bifurcation points can be identified as 
exceptional points: traversing circles around them in the complex-extended parameter 
plane, the eigenvalues are permuted, which is a clear signature of exceptional points. 



4. NONLINEAR DYNAMICS OF BOSE-EINSTEIN 
CONDENSATES WITH ATOMIC LONG-RANGE 
INTERACTIONS 



Starting point of accurate numerical calculations are the time-dependent Gross- 
Pitaevskii equations ([T]) and Q. For variational calculations one makes use of the fact 
that these equations follow from the variational principle ||z^(0 — H'^f{t)\\^ = min. 
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FIGURE 1. Bifurcations of the particle-number scaled chemical potential. Left: 1/r interaction, for 
different trap frequencies 7 (in units of A^^); solid curves: accurate numerical calculation, dashed curves: 
variational calculation. Right: dipole-dipole interaction, variational results for the geometric mean trap 
frequency N'^y = 3.4 x 10'* used in the experiments of Koch et al. I22j and different values of the trap 
aspect ratio X. 



a/a. 



where the variation is performed with respect to and finally ^ is set equal to Using 
a complex parametrization of the trial wave function ^f{t) = x{X{t)), the variation leads 
to equations of motion for the parameters A (cf. [28]) 
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4.1. Time evolution of condensates with 1/r interaction, variational 

and exact 



For simplicity we consider the case of selftrapping, with no external trap. As one 
can convince oneself, the results can be easily generalized to the case where an ex- 
ternal radially symmetric trap potential is present. We choose a Gaussian trial wave 
function V^(r,?) = exp{i[A{t)r^ + Y{t)]}, where A and 7 are complex functions of time, 
whose equations of motion follow from (|4]). Decomposing A into real and imaginary 
parts, A = Ar + iAi, and replacing them by two other dynamical quantities ll29l l30l . 
q = a/3/(4A,;) = a/ (r^), p = A,-^/3/Aj, converts those equations into the canonical 
equations of motion for p and q that follow from the Hamiltonian 



E=H{q,p) = T + V 
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(5) 



In this way the Gross-Pitaevskii equation is mapped onto the Hamiltonian of a one- 
dimensional classical autonomous system with a nonlinear potential V{q). The potential 
has no extremum for a < = — 3;r/8 ~ —1.18, possesses a saddle point for a = a^T, 
and a maximum and a minimum for a > aa- The critical scattering length corresponds 
to the bifurcation point in the variational calculation. For different values of a (in units 
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FIGURE 2. Phase portraits of the dynamics of the complex width function A(f) associated with the 
Hamihonian (jsjl for attractive 1 jr interaction for different values of the scattering length a, measured in 
units of fl„. Left: a = — 1 > a^'- two stationary states appear as fixed points; middle: a = —1.18 ~ a^'- 
coalescence of the fixed points; right: a = — 1.3 < a^x'- no stationary solutions exist. 



of a„) phase portraits of trajectories moving according to the Hamiltonian ([5]) are shown 
in Fig. |2j 

The linear stability analysis of both the variational and the exact quantum stationary 
solutions proves [|30]| that the state corresponding to the elliptical fixed-point indeed is 
dynamically stable (small perturbations of the state show oscillating behavior), while 
the stationary state corresponding to the hyperbolic fixed-point is dynamically unstable 
(exponential growth of small perturbations). 

This behavior is recovered in exact numerical solutions of the time-dependent Gross- 
Pitaevskii equation with 1 /r interaction, but also new features emerge. The solution is 
carried out [30] using the split-operator technique and fast Fourier transforms. To inves- 
tigate the behavior of condensate wave functions in the vicinity of the exact numerical 
stable and unstable stationary states, we consider condensates which are obtained by de- 
forming the stationary states by i//(r) = / ■ (r ■ /^/^ ) , where / is a numerical stretching 
factor (this choice of the perturbation does not affect the norm of the state). 

In Fig. [3] we show examples of the exact BEC dynamics in the vicinity of the unstable 
and stable stationary states. In Fig. |3]a) we start the time evolution with the numerical 
solution for the unstable stationary state (in the classical picture this corresponds to 
the trajectory starting at the hyperbolic fixed point, see left part of Fig. |2]). Because of 
unavoidable numerical deviations from the theoretically exact unstable state, the wave 
function determined numerically is stationary only for some time but then begins to 
oscillate. Obviously we have started with a state which in the variational picture would 
be located in the elliptical domain close to the hyperbolic fixed point. Note, however, 
that the oscillation is not strictly periodic. By contrast, in Fig. [3] b), where the time 
evolution starts with the unstable stationary state stretched by a factor of / = 1.001, 
as time proceeds the wave function contracts towards the origin, and the condensate 
collapses. In the variational picture this corresponds to a trajectory initially close to 
the hyperbolic fixed point but located on the hyperbolic side. Note that in a realistic 
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FIGURE 3. Time evolution, for attractive 1 /r interaction, of the root-mean-square widths of the con- 
densate wave functions in the vicinity of the unstable (panels a), b), c)) and the stable stationary (panels 
d), e)) state, a): Scaled scattering length (in units of a^) a = -1.0, stretching factor f= 1.00; b): a = -0.85, f 
= -1.001; c): a = -0.85, f = 0.99; d): a = -0.85, f =1.25, and e): a = -0.85, f = 1.01. 



experimental situtation during the collapse further mechanisms have to be taken into 
account, such as inelastic collisions. The inclusion of such mechanisms, however, clearly 
goes beyond the scope of the present paper. 

Figure |3]c) displays a behavior not present in the variational analysis. We start again 
in the vicinity of the unstable stationary state (/ = 0.99) and find that the width of 
the condensate gradually grows and grows. An inspection of the wave function on a 
logarithmic scale shows L30J that wave function amplitude builds up at large distances 
from the origin, giving rise to this behavior. Finally, Fig. [3] d), e) show examples for 
the quantum mechanical time evolution of condensates in the vicinity of the stable 
ground state. For a large stretching factor (panel d)) the condensate is found to oscillate 
and to expand, while for a small stretching factor (panel e)) we find the quasiperiodic 
oscillations that we would expect from the variational analysis. This demonstrates that 
the variational nonlinear dynamics approach is capable of predicting essential features 
of the exact quantum mechanical time behavior of the condensates, but that the quantum 
mechanical behavior is even richer. 



4.2. Dynamics of condensates with dipolar interaction, variational 

We choose a Gaussian trial wave function adapted to the axisymmetric trap geom- 
etry, \j/{p,z,t) = e'(^pf +'^;^ +1'), where the complex width parameters Ap, A^, and the 
complex phase are functions of time. Their dynamical equations follow from the time- 



FIGURE 4. Potential V{qp,qS) in the Hamiltonian dTl) for dipole-dipole long-range interaction. 



dependent variational principle Q. Introducing new variables qp^qz^Pp^Pz "^i^ 

KcAp=pp/{Aqp), lmAp = \/{Aql), Re A, = p,/(4^,), ImA, = 1/(8^2^ (6) 

one finds that their dynamical equations are equivalent to the canonical equations of 
motion belonging to the Hamiltonian 

H = T^V = i + i+^ + 2ylql + -4^ + ^+liq' 



2 2 Iql 'P-P iVl^qjqz H 
1 + ql/ql - 3ql ^ctm ^ q^ / {2qj) - \ / (qj^lqj/q^^ - 4^ 
6^q^pqz[\/ql-2/qf 



Thus the variational ansatz has turned the problem into one corresponding to a two- 
dimensional nonintegrable Hamiltonian system, which will exhibit all the features fa- 
miliar from nonlinear dynamics studies of such systems. From the shape of the potential, 
which is shown in Fig. |4]as a function of the "position" variables qp.qz, these features 
can already be read offqualitatively. At the potential minimum sits the stable station- 
ary ground state (elliptic fixed point), while at the saddle point one finds an unstable 
excited stationary state (hyperbolic fixed point). To quantitatively characterize the dy- 
namics of the variational condensate wave functions we follow the trajectories in the 
four-dimensional configuration space spanned by the coordinates of the real and imag- 
inary parts of Ap and A,. Since the total mean-field energy is a constant of motion the 
trajectories are confined to three-dimensional hyperplanes, and their behavior can most 
conveniently be visualized by two-dimensional Poincare surfaces of section defined by 
requiring one of the coordinates to assume a fixed value. 

We consider Poincare surfaces of section defined by the condition that the imaginary 
part of A^{t) is zero. Each time the trajectory crosses the plane Im = 0, the real and 
imaginary parts of Ap{t) = Ap(?) + iA'p{t) are recorded. In Fig. [s] surfaces of section 
are plotted for five different, increasing, values of the mean-field energy. The physical 
parameters of the experiment of Koch et al. [|22| are adopted, and the scattering length 
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FIGURE 5. Poincare surfaces of section of the condensate wave functions for dipolar interaction 
represented by their width parameters at the scaled trap frequency N'^j = 3.4 x 10^, aspect ratio A = 6, 
and the scattering length a/od = 0.1. The surfaces of section correspond to increasing values of the mean- 
field energy (in units of ^d): a) NE = 4.5 x 10^ b) NE = 6.00 x 10^ c) and d) NE = 6.24 x 10^ e) 
NE = 9x 10^ f)NE = 6.00 x lO''. 



is fixed to a/ad = 0.1, away from its critical value. At these parameters, the variational 
mean-field energy of the ground state is NEgs = 4.24 x 10^ (in units of Ed) and represents 
the local minimum on the two-dimensional mean-field energy landscape, plotted as 
a function of the width parameters. The variational energy of the second, unstable, 
stationary state at these experimental parameters is A^^es = 6.24 x 10^, it corresponds 
to the saddle point on the mean-field energy surface. Between these two energy values 
the motion on the trajectories is bound, while for energies above the saddle-point energy 
the motion on the trajectories can become unbound: once the saddle point is traversed 
by a trajectory Ap{t), A^{t), the parameters run to infinity, meaning a shrinking of the 
quantum state to vanishing width, i e., a collapse of the condensate takes place. 

The energy in Fig. |5]a) lies slightly above the energy of the stationary ground state. 
The initially stationary state has evolved into a periodic orbit (fixed point in the surface 
of section), corresponding to a state of the condensate whose motion is periodic. The 
oscillations of the width parameters Ap{t) and A^{t) represent oscillatory stretchings of 
the condensate along the p and z directions. The stable periodic orbit in the surface 
of section is surrounded by elliptical, quasi-periodic orbits, representing quasi-periodic 
oscillations of the condensate. As the energy is increased further, in Fig. |5] b), new 
periodic and quasi-periodic orbits are bom, and the motion is still regular. In Fig. |5] 
c) we have reached the saddle-point energy. Now chaotic orbits have appeared in the 
vicinity of the unstable excited stationary state (hyperbolic fixed point). Figure |5] d) 
shows an enlargement of this region in phase space. In contrast to the (quasi-) periodic 
stretching oscillations of the condensate within the elliptical islands, the chaotic motion 
of the parameters describes a condensate which does not yet collapse but whose widths 
fluctuate irregularly. 



In the surfaces of section shown in Fig. [5] e) and f), with mean-field energies well 
above the saddle-point energy, regular islands are still clearly visible. These stable is- 
lands are surrounded by chaotic trajectories. Since ergodic motion along these trajec- 
tories comes close to every point in the configuration space, the chaotic motion sooner 
or later leads to a crossing of the saddle point and then to the collapse of the conden- 
sate wave functions. It can be seen that with growing energy above the saddle point the 
sizes of the stable regions shrink. The kinematically allowed regions surrounding the 
stable islands are hardly recognizable any more since high above the saddle-point en- 
ergy the chaotic motion becomes more and more unbound, and thus trajectories cross 
the Poincare surfaces of section only a few times, if ever, before they escape to infinity 
and collapse takes place. 

It must be stressed, however, that stable islands persist even far above the saddle-point 
energy, implying the existence of quasi-periodically oscillating nondecaying modes of 
dipolar condensate wave functions. 



5. SUMMARY AND CONCLUSIONS 

We have demonstrated that variational forms of the Bose-Einstein condensate wave 
functions convert the condensates via the Gross-Pitaevskii equation into Hamiltonian 
systems that can be studied using the methods of nonlinear dynamics. We have also 
shown that these results serve as a useful guide in the search for nonlinear dynamics 
effects in the numerically accurate quantum calculations of Bose-Einstein condensates 
with long-range interactions. The existence of stable islands as well as chaotic regions 
for excited states of dipolar Bose-Einstein condensates is a result that could be checked 
experimentally. One way of creating the collectively excited states one might think of is 
to prepare the condensate in the ground state, and then to non-adiabatically reduce the 
trap frequencies. 

One might question whether the Gross-Pitaevskii equation is adequate to describe the 
types of complex dynamics discussed in this paper in "real" condensates. For example, 
in the chaotic regime local density maxima might occur for which losses by two-body 
or three -body collisions would have to be taken into account. However, by virtue of 



the scaling laws discussed in Sections 2.1 and 2.2 parameter ranges can always be found 
where the particle densities remain small even in these regimes, and the Gross-Pitaevskii 
equation is applicable. 

The advantage of the simple variational ansatz adopted in this paper is that the analy- 
sis of the nonlinear dynamical properties of Bose-Einstein condensates becomes partic- 
ularly transparent. Numerical quantum calculations to confirm the variational findings 
for dipolar gases and the extensions to structured condensate states OTl l32ll are under 



way. We have already seen in Sec. 4.1[ by comparing variational and accurate numerical 



quantum results for Bose-Einstein condensates with attractive l/r interaction, that the 
nonlinear dynamical properties predicted by the variational calculation were confirmed 
by the full quantum calculations. We therefore have good reason to believe that this will 
also be true once the full quantum calculations of the dynamics of excited condensate 
wave functions of dipolar gases have become available. 
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